﻿#include "worklet/ComputeResidual.h"
#include "worklet/EquationSet.h"

#include "Application.h"
#include "EulerEquation.h"
#include "Mesh.h"
#include "Transient.h"

#include "TopoTables.h"
#include <vtkm/Matrix.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/cont/ArrayHandleGroupVecVariable.h>
#include <vtkm/worklet/WorkletCellNeighborhood.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/WorkletMapTopology.h>

//template<UInt8 NUM_VARIABLES>
//struct DataPack
//{
//  VTKM_EXEC
//  DataPack()
//  {
//    for (size_t i = 0; i < NUM_VARIABLES; i++)
//    {
//      duh[i] = vtkm::TypeTraits<Vec3f>::ZeroInitialization();
//    }
//  }
//  Vec<FloatDefault, NUM_VARIABLES> uh;
//  Vec<Vec3f, NUM_VARIABLES> duh;
//  //Vec3f duh[NUM_VARIABLES];
//  Vec3f inv_term[NUM_VARIABLES]; // 无粘项
//  Vec3f vis_term[NUM_VARIABLES]; // 粘性项
//
//  VTKM_EXEC
//  void InvisTerm()
//  {
//    FloatDefault _gamma = 1.4;
//    Vec3f vel{ uh[1] / uh[0], uh[2] / uh[0], uh[3] / uh[0] };
//
//    FloatDefault p =
//      (_gamma - 1) * (uh[4] - 0.5 * (uh[1] * uh[1] + uh[2] * uh[2] + uh[3] * uh[3]) / uh[0]);
//    FloatDefault h = (uh[4] + p) / uh[0];
//
//    inv_term[0] = uh[0] * vel;
//
//    inv_term[1] = uh[1] * vel;
//    inv_term[1][0] += p;
//
//    inv_term[2] = uh[2] * vel;
//    inv_term[2][1] += p;
//
//    inv_term[3] = uh[3] * vel;
//    inv_term[3][2] += p;
//
//    inv_term[4] = uh[0] * h * vel;
//  }
//
//  VTKM_EXEC
//  void ViscousTerm()
//  {
//    FloatDefault _reynolds = 1000, _prandtl = 0.72;
//    FloatDefault _gamma = 1.4;
//
//    using RealTensor = vtkm::Matrix<FloatDefault, 3, 3>;
//    //using RealTensor = Vec<Vec3f, 3>;
//    auto r = uh[0];
//    auto drho = duh[0];
//    Vec3f vel{ uh[1] / uh[0], uh[2] / uh[0], uh[3] / uh[0] };
//
//
//    RealTensor grad_mom;
//    //{ duh[1], duh[2], duh[3] };
//    grad_mom[0] = duh[1];
//    grad_mom[1] = duh[2];
//    grad_mom[2] = duh[3];
//
//    RealTensor grad_vel;
//    grad_vel[0] = (grad_mom[0] - vel[0] * drho) / r;
//    grad_vel[1] = (grad_mom[1] - vel[1] * drho) / r;
//    grad_vel[2] = (grad_mom[2] - vel[2] * drho) / r;
//
//    RealTensor grad_vel_transpose = vtkm::MatrixTranspose(grad_vel);
//    FloatDefault vel_div = grad_vel(0, 0) + grad_vel(1, 1) + grad_vel(2, 2);
//
//    RealTensor tau;
//    tau[0] = grad_vel[0] + grad_vel_transpose[0];
//    tau[1] = grad_vel[1] + grad_vel_transpose[1];
//    tau[2] = grad_vel[2] + grad_vel_transpose[2];
//
//    tau(0, 0) -= 2. / 3 * vel_div;
//    tau(1, 1) -= 2. / 3 * vel_div;
//    tau(2, 2) -= 2. / 3 * vel_div;
//
//    FloatDefault vis = 1;
//    for (size_t i = 0; i < 3; i++)
//    {
//      for (size_t j = 0; j < 3; j++)
//      {
//        tau(i, j) *= vis / _reynolds;
//      }
//    }
//
//    Vec3f grad_enthalpy =
//      (duh[4] - uh[4] / uh[0] * duh[0]) / r - MatrixMultiply(vtkm::MatrixTranspose(grad_vel), vel);
//
//    for (size_t i = 0; i < 3; i++)
//    {
//      grad_enthalpy[i] *= (vis / _reynolds) * (_gamma / _prandtl);
//    }
//
//    vis_term[0] = vtkm::TypeTraits<Vec3f>::ZeroInitialization();
//    vis_term[1] = tau[0];
//    vis_term[2] = tau[1];
//    vis_term[3] = tau[2];
//    vis_term[4] = MatrixMultiply(tau, vel) + grad_enthalpy;
//  }
//};


namespace Loops
{

struct AdvancedWorklet : public vtkm::worklet::WorkletMapField
{
  VTKM_CONT
  AdvancedWorklet(FloatDefault dt)
    : _dt(dt)
  {
  }

  using ControlSignature = void(FieldInOut var, FieldIn rhs);
  using ExecutionSignature = void(_1, _2);

  template<typename T>
  VTKM_EXEC void operator()(T& u, const T& rhs) const
  {
    u += -_dt * rhs;
  }

  FloatDefault _dt;
};

struct SolutionToArrayWorklet : public vtkm::worklet::WorkletMapField
{

  using ControlSignature = void(FieldIn solution,
                                FieldOut density,
                                FieldOut velocity,
                                FieldOut rhoe);
  using ExecutionSignature = void(_1, _2, _3, _4);

  template<typename SolType>
  VTKM_EXEC void operator()(const SolType& sol, Real& density, Vec3f& momnetum, Real& rhoe) const
  {
    density = sol[0];
    momnetum = Vec3f{ sol[1], sol[2], sol[3] };
    rhoe = sol[4];
  }
};

//struct ComputeResidualWorklet : vtkm::worklet::WorkletCellNeighborhood
//{
//  VTKM_CONT
//  ComputeResidualWorklet(FloatDefault t)
//    : _t(t)
//  {
//  }
//
//  using ControlSignature =
//    void(CellSetIn, FieldInNeighborhood, FieldOut Rhs, FieldIn centroid, FieldIn spacing);
//  using ExecutionSignature = void(Boundary, _2, _3, _4, _5);
//
//  template<typename NeighIn, typename RhsType>
//  VTKM_EXEC void operator()(const vtkm::exec::BoundaryState& boundary,
//                            const NeighIn& value,
//                            RhsType& rhs,
//                            const Vec3f& centroid,
//                            const Vec3f& spacing) const
//  {
//    //rhs = vtkm::TypeTraits<RhsType>::ZeroInitialization();
//    FloatDefault penalty = 0;
//    Vec3f normal{ 0, 0, 0 };
//    DataPack<RhsType::NUM_COMPONENTS> dl, dr, lift;
//    {
//      //auto& ul = value.Get(0, 0, 0);
//      //auto ur =  value.Get(-1, 0, 0);
//      dl.uh = value.Get(0, 0, 0);
//      dr.uh = value.Get(-1, 0, 0);
//      RhsType lflux, rflux;
//      if (boundary.IJK[0] == 0)
//      {
//        Vec3f p = centroid;
//        p[0] = p[0] - 0.5 * spacing[0];
//        //ur = ul;
//        dr.uh = value.Get(0, 0, 0);
//        //dr.uh = valueExactCoutteFlow(_t, p);
//      }
//      dl.InvisTerm();
//      dr.InvisTerm();
//
//      normal = Vec3f{ -1, 0, 0 };
//      for (size_t i = 0; i < RhsType::NUM_COMPONENTS; ++i)
//      {
//        lift.uh[i] = (dl.uh[i] + dr.uh[i]) / 2;
//        lift.duh[i] = (dl.uh[i] - dr.uh[i]) / 2 * normal;
//      }
//
//      RiemannFlux(dl, dr, normal, lflux);
//      lift.ViscousTerm();
//      //std::cout << lift.vis_term[0] << std::endl;
//      //std::cout << lift.vis_term[1] << std::endl;
//      //std::cout << lift.vis_term[2] << std::endl;
//      //std::cout << lift.vis_term[3] << std::endl;
//      //std::cout << lift.vis_term[4] << std::endl;
//
//      for (size_t i = 0; i < RhsType::NUM_COMPONENTS; i++)
//      {
//        lflux[i] += penalty * vtkm::Dot(lift.vis_term[i], normal) / spacing[0];
//      }
//
//
//      dl.uh = value.Get(0, 0, 0);
//      dr.uh = value.Get(1, 0, 0);
//      if (boundary.IJK[0] == boundary.PointDimensions[0] - 1)
//      {
//        Vec3f p = centroid;
//        p[0] = p[0] + 0.5 * spacing[0];
//        //ur = ul;
//        dr.uh = value.Get(0, 0, 0);
//        //dr.uh = valueExactCoutteFlow(_t, p);
//      }
//      dl.InvisTerm();
//      dr.InvisTerm();
//
//
//      normal = Vec3f{ 1, 0, 0 };
//      for (size_t i = 0; i < RhsType::NUM_COMPONENTS; ++i)
//      {
//        lift.uh[i] = (dl.uh[i] + dr.uh[i]) / 2;
//        lift.duh[i] = (dl.uh[i] - dr.uh[i]) / 2 * normal;
//      }
//
//      RiemannFlux(dl, dr, normal, rflux);
//      lift.ViscousTerm();
//      for (size_t i = 0; i < RhsType::NUM_COMPONENTS; i++)
//      {
//        rflux[i] += penalty * vtkm::Dot(lift.vis_term[i], normal) / spacing[0];
//      }
//
//      //RiemannFlux(ul, ur, Vec3f{ 1, 0, 0 }, rflux);
//      rhs = (rflux + lflux) / spacing[0];
//    }
//
//    {
//      //auto& ul = value.Get(0, 0, 0);
//      //auto ur = value.Get(0, -1, 0);
//      dl.uh = value.Get(0, 0, 0);
//      dr.uh = value.Get(0, -1, 0);
//      RhsType lflux, rflux;
//      if (boundary.IJK[1] == 0)
//      {
//        Vec3f p = centroid;
//        p[1] = p[1] - 0.5 * spacing[1];
//        dr.uh = value.Get(0, 0, 0);
//        //dr.uh = valueExactCoutteFlow(_t, p);
//      }
//
//      //dl.uh = ul;
//      //dr.uh = ur;
//      dl.InvisTerm();
//      dr.InvisTerm();
//
//      normal = Vec3f{ 0, -1, 0 };
//      for (size_t i = 0; i < RhsType::NUM_COMPONENTS; ++i)
//      {
//        lift.uh[i] = (dl.uh[i] + dr.uh[i]) / 2;
//        lift.duh[i] = (dl.uh[i] - dr.uh[i]) / 2 * normal;
//      }
//
//      RiemannFlux(dl, dr, normal, lflux);
//      lift.ViscousTerm();
//      for (size_t i = 0; i < RhsType::NUM_COMPONENTS; i++)
//      {
//        lflux[i] += penalty * vtkm::Dot(lift.vis_term[i], normal) / spacing[0];
//      }
//
//      //ur = value.Get(0, 1, 0);
//      dl.uh = value.Get(0, 0, 0);
//      dr.uh = value.Get(0, 1, 0);
//      if (boundary.IJK[1] == boundary.PointDimensions[1] - 1)
//      {
//        Vec3f p = centroid;
//        p[1] = p[1] + 0.5 * spacing[1];
//        dr.uh = value.Get(0, 0, 0);
//        //dr.uh = valueExactCoutteFlow(_t, p);
//      }
//      dl.InvisTerm();
//      dr.InvisTerm();
//
//      normal = Vec3f{ 0, 1, 0 };
//      for (size_t i = 0; i < RhsType::NUM_COMPONENTS; ++i)
//      {
//        lift.uh[i] = (dl.uh[i] + dr.uh[i]) / 2;
//        lift.duh[i] = (dl.uh[i] - dr.uh[i]) / 2 * normal;
//      }
//
//      RiemannFlux(dl, dr, normal, rflux);
//      lift.ViscousTerm();
//      for (size_t i = 0; i < RhsType::NUM_COMPONENTS; i++)
//      {
//        rflux[i] += penalty * vtkm::Dot(lift.vis_term[i], normal) / spacing[0];
//      }
//
//      rhs += (rflux + lflux) / spacing[1];
//    }
//  }
//
//  VTKM_EXEC
//  vtkm::Vec<FloatDefault, 5> valueExactIsoVortex(FloatDefault t, const Vec3f& p) const
//  {
//    Vec3f _vel{ 1, 1, 0 };
//    FloatDefault x = p[0] - _vel[0] * t;
//    FloatDefault y = p[1] - _vel[1] * t;
//    FloatDefault z = p[2] - _vel[2] * t;
//    FloatDefault gam = 1.4, gamm1 = gam - 1, epi = 5.;
//    FloatDefault xb, yb, r2;
//
//    FloatDefault PI = vtkm::Pi();
//    xb = x + 5;
//    yb = y + 5;
//    r2 = xb * xb + yb * yb;
//
//    FloatDefault u = _vel[0] + epi / (2. * PI) * exp(0.5 * (1. - r2)) * (-yb);
//    FloatDefault v = _vel[1] + epi / (2. * PI) * exp(0.5 * (1. - r2)) * xb;
//    FloatDefault T = 1. - gamm1 * epi * epi / (8 * gam * PI * PI) * exp(1. - r2);
//    FloatDefault density = pow(T, 1 / gamm1);
//    FloatDefault pre = pow(density, gam);
//
//    vtkm::Vec<FloatDefault, 5> value;
//    value[0] = density;
//    value[1] = density * u;
//    value[2] = density * v;
//    value[3] = density * 0;
//    value[4] = pre / gamm1 + 0.5 * density * (u * u + v * v);
//
//    return value;
//  }
//
//  VTKM_EXEC
//  vtkm::Vec<FloatDefault, 5> valueExactCoutteFlow(FloatDefault t, const Vec3f& p) const
//  {
//    FloatDefault _gamma = 1.4, _prandtl = 0.72, _mach = 0.2;
//
//    FloatDefault x = p[0];
//    FloatDefault y = p[1];
//    FloatDefault z = p[2];
//
//    FloatDefault T2 = 0.85, T1 = 0.8;
//    FloatDefault tem =
//      T1 + (T2 - T1) * y / 2 + 0.5 * _prandtl * (_gamma - 1) * _mach * _mach * y / 2 * (1 - y / 2);
//    FloatDefault pre = 1. / (_gamma * _mach * _mach);
//    FloatDefault density = pre * _gamma * _mach * _mach / tem;
//    FloatDefault u = y / 2.0, v = 0, w = 0;
//    Vec3f mom{ density * u, density * v, density * w };
//    FloatDefault rhoe = pre / (_gamma - 1) + 0.5 * density * (u * u + v * v + w * w);
//
//    vtkm::Vec<FloatDefault, 5> value;
//    value[0] = density;
//    value[1] = density * u;
//    value[2] = density * v;
//    value[3] = density * w;
//    value[4] = rhoe;
//
//    return value;
//  }
//  template<typename DataPackType, typename VType>
//  VTKM_EXEC void RiemannFlux(const DataPackType& dl,
//                             const DataPackType& dr,
//                             const Vec3f& normal,
//                             VType& riemann_flux) const
//  {
//    for (IdComponent eq = 0; eq < VType::NUM_COMPONENTS; eq++)
//    {
//      riemann_flux[eq] =
//        vtkm::Dot(dl.inv_term[eq] + dr.inv_term[eq], normal) / 2 + 0.6 * (dl.uh[eq] - dr.uh[eq]);
//    }
//  }
//  FloatDefault _t;
//};

//struct ComputeResidual_ForwardStep : vtkm::worklet::WorkletCellNeighborhood
//{
//  VTKM_CONT
//  ComputeResidual_ForwardStep(FloatDefault t)
//    : _t(t)
//  {
//  }
//
//  using ControlSignature = void(CellSetIn,
//                                FieldInNeighborhood solution,
//                                FieldOut Rhs,
//                                FieldIn centroid,
//                                FieldIn spacing,
//                                FieldIn blk_id,
//                                FieldIn bnd_id);
//  using ExecutionSignature = void(Boundary, _2, _3, _4, _5, _6, _7);
//
//  template<typename NeighIn, typename RhsType>
//  VTKM_EXEC void operator()(const vtkm::exec::BoundaryState& boundary,
//                            const NeighIn& value,
//                            RhsType& rhs,
//                            const Vec3f& centroid,
//                            const Vec3f& spacing,
//                            const vtkm::Int16& blk_id,
//                            const Vec<vtkm::Int16, 4>& bnd_id) const
//  {
//    if (blk_id == -1)
//    {
//      rhs[0] = 0;
//      rhs[1] = 0;
//      rhs[2] = 0;
//      rhs[3] = 0;
//      rhs[4] = 0;
//      return;
//    }
//    Vec3f normal{ 0, 0, 0 };
//    Vec3f normal_wall{ 0, 0, 0 };
//    DataPack<RhsType::NUM_COMPONENTS> dl, dr;
//    dl.uh = value.Get(0, 0, 0);
//    dl.InvisTerm();
//    {
//      normal = Vec3f{ -1, 0, 0 };
//      dr.uh = value.Get(-1, 0, 0);
//      RhsType lflux, rflux;
//      Vec3f p = centroid;
//      p[0] = p[0] - 0.5 * spacing[0];
//      normal_wall = p / vtkm::Magnitude(p);
//      BoundaryCondition(bnd_id[3], _t, p, normal_wall, dl, dr);
//      dr.InvisTerm();
//      RiemannFlux(dl, dr, normal, lflux);
//
//      dr.uh = value.Get(1, 0, 0);
//      normal = Vec3f{ 1, 0, 0 };
//      p = centroid;
//      p[0] = p[0] + 0.5 * spacing[0];
//      normal_wall = p / vtkm::Magnitude(p);
//      BoundaryCondition(bnd_id[1], _t, p, normal, dl, dr);
//      dr.InvisTerm();
//      RiemannFlux(dl, dr, normal, rflux);
//
//      rhs = (rflux + lflux) / spacing[0];
//    }
//
//    {
//      dr.uh = value.Get(0, -1, 0);
//      RhsType lflux, rflux;
//      normal = Vec3f{ 0, -1, 0 };
//      Vec3f p = centroid;
//      p[1] = p[1] - 0.5 * spacing[1];
//      normal_wall = p / vtkm::Magnitude(p);
//      BoundaryCondition(bnd_id[0], _t, p, normal, dl, dr);
//      dr.InvisTerm();
//      RiemannFlux(dl, dr, normal, lflux);
//
//      dr.uh = value.Get(0, 1, 0);
//      normal = Vec3f{ 0, 1, 0 };
//      p = centroid;
//      p[1] = p[1] + 0.5 * spacing[1];
//      normal_wall = p / vtkm::Magnitude(p);
//      BoundaryCondition(bnd_id[2], _t, p, normal, dl, dr);
//      dr.InvisTerm();
//      RiemannFlux(dl, dr, normal, rflux);
//
//      rhs += (rflux + lflux) / spacing[1];
//    }
//  }
//
//  template<typename DataPackType, typename VType>
//  VTKM_EXEC void RiemannFlux(const DataPackType& dl,
//                             const DataPackType& dr,
//                             const Vec3f& normal,
//                             VType& riemann_flux) const
//  {
//    for (IdComponent eq = 0; eq < VType::NUM_COMPONENTS; eq++)
//    {
//      riemann_flux[eq] =
//        vtkm::Dot(dl.inv_term[eq] + dr.inv_term[eq], normal) / 2 + 1.5 * (dl.uh[eq] - dr.uh[eq]);
//    }
//  }
//
//  template<typename DataPackType>
//  VTKM_EXEC void BoundaryCondition(const vtkm::Int16& bid,
//                                   const FloatDefault& t,
//                                   const Vec3f& p,
//                                   const Vec3f& normal,
//                                   const DataPackType& dl,
//                                   DataPackType& dr) const
//  {
//    if (bid == 1) // 出口
//    {
//      auto& ul = dl.uh;
//      auto& ur = dr.uh;
//      ur[0] = ul[0];
//      ur[1] = ul[1];
//      ur[2] = ul[2];
//      ur[3] = ul[3];
//      ur[4] = ul[4];
//    }
//    if (bid == 3) // 进口
//    {
//      using ValueType = Vec<FloatDefault, 5>;
//      ValueType primary_value;
//      primary_value = ValueType{ 1.4, 3, 0, 0, 1 };
//
//      FloatDefault gam = 1.4, gamm1 = gam - 1;
//      ValueType value;
//      auto r = primary_value[0];
//      auto u = primary_value[1];
//      auto v = primary_value[2];
//      auto w = primary_value[3];
//      auto pre = primary_value[4];
//
//      auto& ur = dr.uh;
//
//      ur[0] = r;
//      ur[1] = r * u;
//      ur[2] = r * v;
//      ur[3] = r * w;
//      ur[4] = pre / gamm1 + 0.5 * r * (u * u + v * v);
//    }
//    if (bid == 0 || bid == 2 || bid == 4) //wall
//    {
//      auto& ul = dl.uh;
//      auto& ur = dr.uh;
//      Vec3f mom{ ul[1], ul[2], ul[3] };
//      Real vn = vtkm::Dot(mom, normal);
//      ur[0] = ul[0];
//      ur[1] = ul[1] - 2 * vn * normal[0];
//      ur[2] = ul[2] - 2 * vn * normal[1];
//      ur[3] = ul[3] - 2 * vn * normal[2];
//      ur[4] = ul[4];
//    }
//  }
//
//  FloatDefault _t;
//};


//struct ComputeResidual_IsoVortex : vtkm::worklet::WorkletCellNeighborhood
//{
//  VTKM_CONT
//  ComputeResidual_IsoVortex(FloatDefault t)
//    : _t(t)
//  {
//  }
//
//  using ControlSignature =
//    void(CellSetIn, FieldInNeighborhood, FieldOut Rhs, FieldIn centroid, FieldIn spacing);
//  using ExecutionSignature = void(Boundary, _2, _3, _4, _5);
//
//  template<typename NeighIn, typename RhsType>
//  VTKM_EXEC void operator()(const vtkm::exec::BoundaryState& boundary,
//                            const NeighIn& value,
//                            RhsType& rhs,
//                            const Vec3f& centroid,
//                            const Vec3f& spacing) const
//  {
//    Vec3f normal{ 0, 0, 0 };
//    DataPack<RhsType::NUM_COMPONENTS> dl, dr;
//    Real k = 1. / 3;
//    using ValueType = Vec<FloatDefault, 5>;
//
//
//    { // x方向
//      normal = Vec3f{ -1, 0, 0 };
//      //dl.uh = InterLeft(value.Get(-1, 0, 0), value.Get(0, 0, 0), value.Get(1, 0, 0));
//      //dr.uh = InterRight(value.Get(-2, 0, 0), value.Get(-1, 0, 0), value.Get(0, 0, 0));
//      dl.uh = value.Get(0, 0, 0);
//      dr.uh = value.Get(-1, 0, 0);
//      RhsType lflux, rflux;
//      Vec3f p = centroid;
//      p[0] = p[0] - 0.5 * spacing[0];
//      //BoundaryCondition(3, _t, p, normal, dl, dr);
//      dl.InvisTerm();
//      dr.InvisTerm();
//      RiemannFlux(dl, dr, normal, lflux);
//
//      //dl.uh = InterRight(value.Get(-1, 0, 0), value.Get(0, 0, 0), value.Get(1, 0, 0));
//      //dr.uh = InterLeft(value.Get(0, 0, 0), value.Get(1, 0, 0), value.Get(2, 0, 0));
//      dl.uh = value.Get(0, 0, 0);
//      dr.uh = value.Get(1, 0, 0);
//      normal = Vec3f{ 1, 0, 0 };
//      p = centroid;
//      p[0] = p[0] + 0.5 * spacing[0];
//      //BoundaryCondition(1, _t, p, normal, dl, dr);
//      dl.InvisTerm();
//      dr.InvisTerm();
//      RiemannFlux(dl, dr, normal, rflux);
//
//      rhs = (rflux + lflux) / spacing[0];
//    }
//    //std::cout << rhs << std::endl;
//
//    { // y方向
//      normal = Vec3f{ 0, -1, 0 };
//      //dl.uh = InterLeft(value.Get(0, -1, 0), value.Get(0, 0, 0), value.Get(0, 1, 0));
//      //dr.uh = InterRight(value.Get(0, -2, 0), value.Get(0, -1, 0), value.Get(0, 0, 0));
//      dl.uh = value.Get(0, 0, 0);
//      dr.uh = value.Get(0, -1, 0);
//      RhsType lflux, rflux;
//      Vec3f p = centroid;
//      p[1] = p[1] - 0.5 * spacing[1];
//      //BoundaryCondition(0, _t, p, normal, dl, dr);
//      dl.InvisTerm();
//      dr.InvisTerm();
//      RiemannFlux(dl, dr, normal, lflux);
//
//      normal = Vec3f{ 0, 1, 0 };
//      //dl.uh = InterRight(value.Get(0, -1, 0), value.Get(0, 0, 0), value.Get(0, 1, 0));
//      //dr.uh = InterLeft(value.Get(0, 0, 0), value.Get(0, 1, 0), value.Get(0, 2, 0));
//      dl.uh = value.Get(0, 0, 0);
//      dr.uh = value.Get(0, 1, 0);
//      p = centroid;
//      p[1] = p[1] + 0.5 * spacing[1];
//      //BoundaryCondition(2, _t, p, normal, dl, dr);
//      dl.InvisTerm();
//      dr.InvisTerm();
//      RiemannFlux(dl, dr, normal, rflux);
//
//      rhs += (rflux + lflux) / spacing[1];
//    }
//
//    //{ // z方向
//    //  dl.uh = value.Get(0, 0, 0);
//    //  dl.InvisTerm();
//    //  //MUSCLLeft(dl, value.Get(0, 0, 0), value.Get(0, 0, 0), value.Get(0, 0, 0));
//    //  dr.uh = value.Get(0, 0, -1);
//    //  RhsType lflux, rflux;
//    //  normal = Vec3f{ 0, 0, -1 };
//    //  Vec3f p = centroid;
//    //  p[2] = p[2] - 0.5 * spacing[2];
//    //  BoundaryCondition(0, _t, p, normal, dl, dr);
//    //  dr.InvisTerm();
//    //  RiemannFlux(dl, dr, normal, lflux);
//
//    //  dl.uh = value.Get(0, 0, 0);
//    //  dl.InvisTerm();
//    //  dr.uh = value.Get(0, 0, 1);
//    //  normal = Vec3f{ 0, 0, 1 };
//    //  p = centroid;
//    //  p[2] = p[2] + 0.5 * spacing[2];
//    //  BoundaryCondition(2, _t, p, normal, dl, dr);
//    //  dr.InvisTerm();
//    //  RiemannFlux(dl, dr, normal, rflux);
//
//    //  rhs += (rflux + lflux) / spacing[2];
//    //}
//  }
//
//  template<typename VType>
//  VTKM_EXEC VType FirstOrderLeft(VType& _l, VType& _0, VType& _r) const
//  {
//    return _0;
//  }
//
//  template<typename VType>
//  VTKM_EXEC VType FirstOrderRight(VType& _l, VType& _0, VType& _r) const
//  {
//    return _0;
//  }
//
//  template<typename VType>
//  VTKM_EXEC VType SecondOrderLeft(VType& _l, VType& _0, VType& _r) const
//  {
//    return 2 * _0 - _r;
//  }
//
//  template<typename VType>
//  VTKM_EXEC VType SecondOrderRight(VType& _l, VType& _0, VType& _r) const
//  {
//    return 2 * _0 - _l;
//  }
//
//  template<typename VType>
//  VTKM_EXEC VType MUSCLLeft(VType& _l, VType& _0, VType& _r) const
//  {
//    auto delta_l = _0 - _l;
//    auto delta_r = _r - _0;
//    Real k = 1 / .3;
//    return _0 - 1. / 4 * ((1 - k) * delta_l + (1 + k) * delta_r);
//  }
//
//  template<typename VType>
//  VTKM_EXEC VType MUSCLRight(VType& _l, VType& _0, VType& _r) const
//  {
//    auto delta_l = _0 - _l;
//    auto delta_r = _r - _0;
//    Real k = 1 / .3;
//    return _0 + 1. / 4 * ((1 - k) * delta_l + (1 + k) * delta_r);
//  }
//
//  template<typename VType>
//  VTKM_EXEC VType InterLeft(VType& _l, VType& _0, VType& _r) const
//  {
//    return FirstOrderLeft(_l, _0, _r);
//    //return SecondOrderLeft(_l, _0, _r);
//    //return MUSCLLeft(_l, _0, _r);
//  }
//
//  template<typename VType>
//  VTKM_EXEC VType InterRight(VType& _l, VType& _0, VType& _r) const
//  {
//    return FirstOrderRight(_l, _0, _r);
//    //return SecondOrderRight(_l, _0, _r);
//    //return MUSCLRight(_l, _0, _r);
//  }
//
//  template<typename DataPackType, typename VType>
//  VTKM_EXEC void RiemannFlux(const DataPackType& dl,
//                             const DataPackType& dr,
//                             const Vec3f& normal,
//                             VType& riemann_flux) const
//  {
//    for (IdComponent eq = 0; eq < VType::NUM_COMPONENTS; eq++)
//    {
//      riemann_flux[eq] =
//        vtkm::Dot(dl.inv_term[eq] + dr.inv_term[eq], normal) / 2 + 2 * (dl.uh[eq] - dr.uh[eq]);
//    }
//  }
//
//  template<typename DataPackType>
//  VTKM_EXEC void BoundaryCondition(const vtkm::Int16& bid,
//                                   const FloatDefault& t,
//                                   const Vec3f& p,
//                                   const Vec3f& normal,
//                                   const DataPackType& dl,
//                                   DataPackType& dr) const
//  {
//    //dr.uh = IsoVortexExact(t, p);
//    dr.uh = dl.uh;
//  }
//
//  using ValueType = Vec<FloatDefault, 5>;
//  VTKM_EXEC ValueType IsoVortexExact(const FloatDefault t, const Vec3f& p) const
//  {
//    ValueType value;
//    Vec3f _vel{ 1, 1, 0 };
//    FloatDefault x = p[0] - _vel[0] * t;
//    FloatDefault y = p[1] - _vel[1] * t;
//    FloatDefault z = p[2] - _vel[2] * t;
//    FloatDefault gam = 1.4, gamm1 = gam - 1, epi = 5.;
//    FloatDefault xb, yb, r2;
//
//    FloatDefault PI = vtkm::Pi();
//    xb = x + 5;
//    yb = y + 5;
//    r2 = xb * xb + yb * yb;
//
//    FloatDefault u = _vel[0] + epi / (2. * PI) * exp(0.5 * (1. - r2)) * (-yb);
//    FloatDefault v = _vel[1] + epi / (2. * PI) * exp(0.5 * (1. - r2)) * xb;
//    FloatDefault T = 1. - gamm1 * epi * epi / (8 * gam * PI * PI) * exp(1. - r2);
//    FloatDefault density = pow(T, 1 / gamm1);
//    FloatDefault pre = pow(density, gam);
//
//    value[0] = density;
//    value[1] = density * u;
//    value[2] = density * v;
//    value[3] = density * 0;
//    value[4] = pre / gamm1 + 0.5 * density * (u * u + v * v);
//
//    return value;
//  }
//  FloatDefault _t;
//};



struct ComputeResidualWorkletPassFlow : vtkm::worklet::WorkletVisitCellsWithPoints
{
  VTKM_CONT
  ComputeResidualWorkletPassFlow(FloatDefault t)
    : _t(t)
  {
  }

  using ControlSignature = void(CellSetIn cell_set_in,
                                FieldInPoint coords,
                                FieldInCell cell_link,
                                WholeArrayIn sln,
                                FieldOutCell rhs,
                                ExecObject es);

  using ExecutionSignature = void(CellShape, InputIndex, _2, _3, _4, _5, _6);

  template<typename PointCoordVecType>
  VTKM_EXEC Vec3f ComputeDirectedArea(const vtkm::CellShapeTagGeneric& shape,
                                      const PointCoordVecType& points) const // 有向面积
  {
    Vec3f directed_area{ 1, 1, 1 };
    switch (shape.Id)
    {
      case vtkm::CELL_SHAPE_VERTEX:
      {
        this->RaiseError("1D vertex normal unsolved.");
      }
      break;
      case vtkm::CELL_SHAPE_LINE:
      {
        using CellShape = vtkm::CellShapeIdToTag<vtkm::CELL_SHAPE_LINE>::Tag;
        auto _u_axis = points[1] - points[0];
        directed_area = { _u_axis[1], -_u_axis[0], 0 };
        //directed_area = vtkm::Cross(points[1] - points[0], Vec3f{ 0,0,1 });
      }
      break;

      case vtkm::CELL_SHAPE_TRIANGLE:
      {
        directed_area = vtkm::Cross(points[0] - points[1], points[2] - points[1]);
        directed_area = directed_area / 2;
      }
      break;
      case vtkm::CELL_SHAPE_QUAD:
      {
        directed_area = vtkm::Cross(points[0] - points[1], points[2] - points[1]);
        directed_area += vtkm::Cross(points[2] - points[3], points[0] - points[3]);
        directed_area = directed_area / 2;
      }
      break;
      default:
        this->RaiseError("unknown cell type %d");
        break;
    }
    return directed_area;
  }

  template<typename PointCoordVecType>
  VTKM_EXEC Vec3f FaceCenter(const vtkm::CellShapeTagGeneric& shape,
                             const PointCoordVecType& pts) const // 有向面积
  {
    Vec3f centroid{ 1, 1, 1 };
    switch (shape.Id)
    {
      case vtkm::CELL_SHAPE_VERTEX:
      {
        this->RaiseError("1D vertex normal unsolved.");
      }
      break;
      case vtkm::CELL_SHAPE_LINE:
      {
        centroid = (pts[0] + pts[1]) / 2;
      }
      break;

      case vtkm::CELL_SHAPE_TRIANGLE:
      {
        centroid = (pts[0] + pts[1] + pts[2]) / 3;
      }
      break;
      case vtkm::CELL_SHAPE_QUAD:
      {
        centroid = (pts[0] + pts[1] + pts[2] + pts[3]) / 4;
      }
      break;
      default:
        this->RaiseError("unknown cell type %d");
        break;
    }
    return centroid;
  }


  template<typename CellShapeTag,
           typename CoordType,
           typename CellLink,
           typename SlnType,
           typename RhsType/*,
           typename EquationType*/>
  VTKM_EXEC void operator()(CellShapeTag& cell_shape,
                            Id& cell_id,
                            const CoordType& pts,
                            const CellLink& cell_link,
                            const SlnType& sln,
                            RhsType& rhs,
                            EquationSet& es) const
  {
    auto num_faces = TopoTables{}.NumSides(cell_shape.Id);
    Real volume = 0;
    rhs = vtkm::TypeTraits<RhsType>::ZeroInitialization();
    using ValueType = RhsType;
    ValueType ul, ur;
    RhsType flux;

    Vec3f f2n[4];
    for (IdComponent side = 0; side < num_faces; side++)
    {
      UInt8 face_shape;
      CellFaceShape(side, cell_shape, face_shape);
      IdComponent n_f2n;
      CellFaceNumberOfPoints(side, cell_shape, n_f2n);

      for (IdComponent i = 0; i < n_f2n; i++)
      {
        IdComponent result;
        CellFaceLocalIndex(i, side, cell_shape, result);
        f2n[i] = pts[result];
      }

      ValueType neighbor_value;
      Id neighbor_cell_id = cell_link[side].Get().first;
      Vec3f centroid = FaceCenter(face_shape, f2n);

      Vec3f normal = ComputeDirectedArea(face_shape, f2n);
      const Real face_area = vtkm::Magnitude(normal);
      normal = -normal / face_area;

      ul = sln.Get(cell_id);
      ur = sln.Get(cell_id);
      if (neighbor_cell_id < 0) // 边界条件
      {
        auto bid = cell_link[side].Get().second;
        ur = es.BoundaryCondition(bid, ul, normal, centroid);
      }
      else
      {
        ur = sln.Get(neighbor_cell_id);
      }
      es.RiemannFlux(ul, ur, normal, flux);
      rhs += face_area * flux;
      volume += normal[1] * centroid[1] * face_area;
    }
    ////std::cout << volume << std::endl;
    rhs = rhs / (volume);
  }
  Real _t;
};


struct VolumeWorklet : vtkm::worklet::WorkletVisitCellsWithPoints
{

  using ControlSignature = void(CellSetIn cell_set_in,
                                FieldInPoint coords,
                                FieldInCell cell_link,
                                FieldOutCell volume);

  using ExecutionSignature = void(CellShape, _2, _3, _4);

  template<typename PointCoordVecType>
  VTKM_EXEC Vec3f ComputeDirectedArea(const vtkm::CellShapeTagGeneric& shape,
                                      const PointCoordVecType& points) const // 有向面积
  {
    Vec3f directed_area{ 1, 1, 1 };
    switch (shape.Id)
    {
      case vtkm::CELL_SHAPE_VERTEX:
      {
        this->RaiseError("1D vertex normal unsolved.");
      }
      break;
      case vtkm::CELL_SHAPE_LINE:
      {
        using CellShape = vtkm::CellShapeIdToTag<vtkm::CELL_SHAPE_LINE>::Tag;
        auto _u_axis = points[1] - points[0];
        directed_area = { _u_axis[1], -_u_axis[0], 0 };
        //directed_area = vtkm::Cross(points[1] - points[0], Vec3f{ 0,0,1 });
      }
      break;

      case vtkm::CELL_SHAPE_TRIANGLE:
      {
        directed_area = vtkm::Cross(points[0] - points[1], points[2] - points[1]);
        directed_area = directed_area / 2;
      }
      break;
      case vtkm::CELL_SHAPE_QUAD:
      {
        directed_area = vtkm::Cross(points[0] - points[1], points[2] - points[1]);
      }
      break;
      default:
        this->RaiseError("unknown cell type %d");
        break;
    }
    return directed_area;
  }
  template<typename PointCoordVecType>
  VTKM_EXEC Vec3f FaceCenter(const vtkm::CellShapeTagGeneric& shape,
                             const PointCoordVecType& pts) const // 有向面积
  {
    Vec3f centroid{ 1, 1, 1 };
    switch (shape.Id)
    {
      case vtkm::CELL_SHAPE_VERTEX:
      {
        this->RaiseError("1D vertex normal unsolved.");
      }
      break;
      case vtkm::CELL_SHAPE_LINE:
      {
        centroid = (pts[0] + pts[1]) / 2;
      }
      break;

      case vtkm::CELL_SHAPE_TRIANGLE:
      {
        centroid = (pts[0] + pts[1] + pts[2]) / 3;
      }
      break;
      case vtkm::CELL_SHAPE_QUAD:
      {
        centroid = (pts[0] + pts[1] + pts[2] + pts[3]) / 4;
      }
      break;
      default:
        this->RaiseError("unknown cell type %d");
        break;
    }
    return centroid;
  }

  template<typename CellShapeTag, typename CoordType, typename CellLink>
  VTKM_EXEC void operator()(CellShapeTag cell_shape,
                            const CoordType& pts,
                            const CellLink& cell_link,
                            Real& volume) const
  {
    auto num_faces = TopoTables{}.NumSides(cell_shape.Id);
    volume = 0;
    Vec3f f2n[4];
    for (IdComponent side = 0; side < num_faces; side++)
    {
      UInt8 face_shape;
      CellFaceShape(side, cell_shape, face_shape);
      IdComponent n_f2n;
      CellFaceNumberOfPoints(side, cell_shape, n_f2n);
      for (IdComponent i = 0; i < n_f2n; i++)
      {
        IdComponent result;
        CellFaceLocalIndex(i, side, cell_shape, result);
        f2n[i] = pts[result];
      }

      Vec3f normal = ComputeDirectedArea(face_shape, f2n);
      const Real face_area = vtkm::Magnitude(normal);
      //std::cout << face_area << std::endl;
      normal = normal / face_area;
      Vec3f centroid = FaceCenter(face_shape, f2n);
      volume += normal[0] * centroid[0] * face_area;
      //if (side == 0)
      //volume = face_area;
    }
    std::cout << volume << std::endl;
  }
};

void ExplictFirtOrder(EulerEquation& eq)
{
  auto& mesh = eq.GetMesh();
  auto& sln = eq.GetSolution();
  auto& rhs = eq.GetRhs();
  auto& transient = static_cast<Transient&>(*(eq.App().GetExecutioner()));
  Real dt = transient.Dt();

  Invoker{}(AdvancedWorklet{ dt }, sln, rhs);
}

void SolutionToArray(EulerEquation& eq)
{
  auto& mesh = eq.GetMesh();
  auto& sln = eq.GetSolution();
  auto& rhs = eq.GetRhs();
  auto density = eq.GetVariableArray<Real>("density");
  auto momentum = eq.GetVariableArray<Vec3f>("momentum");
  auto rhoe = eq.GetVariableArray<Real>("rhoe");
  Invoker{}(SolutionToArrayWorklet{}, sln, density, momentum, rhoe);
}

void ComputeResidual(EulerEquation& sys)
{
  auto& transient = static_cast<Transient&>(*(sys.App().GetExecutioner()));
  Real current_time = transient.CurrentTime();

  Timer timer;
  timer.Start();
  auto& mesh = sys.GetMesh();
  vtkm::cont::Invoker{}(ComputeResidualWorkletPassFlow{ current_time },
                        mesh.GetCellSet(),
                        mesh.GetCoordinate(),
                        mesh.GetCellLinks(),
                        sys.GetSolution(),
                        sys.GetRhs(),
                        sys.GetEquationSet());
}

void ComputeResidualIsoVortex(EulerEquation& eq)
{
  auto& transient = static_cast<Transient&>(*(eq.App().GetExecutioner()));
  Real current_time = transient.CurrentTime();

  Timer timer;
  timer.Start();
  const auto& centroid = eq.GetMesh().GetFieldArray<Vec3f>("center");
  const auto& spacing = eq.GetMesh().GetFieldArray<Vec3f>("spacing");
  auto& sln = eq.GetSolution();
  auto& rhs = eq.GetRhs();

  //vtkm::cont::Invoker{}(ComputeResidual_IsoVortex{ current_time },
  //                      eq.GetMesh().GetCellSet(),
  //                      sln,
  //                      rhs,
  //                      centroid,
  //                      spacing);
}

void ComputeResidualIsoVortexOnUnstructuredMesh(EulerEquation& eq)
{
  //auto& transient = static_cast<Transient&>(*(eq.App().GetExecutioner()));
  //Real current_time = transient.CurrentTime();

  //Timer timer;
  //timer.Start();
  //auto& sln = eq.GetSolution();
  //auto& rhs = eq.GetRhs();
  //auto& mesh = eq.GetMesh();

  //auto& face_num_offset = mesh.GetFaceNumOffset();
  //auto& cell_links =
  //  vtkm::cont::make_ArrayHandleGroupVecVariable(mesh.GetCellLinks(), face_num_offset);
  //auto& coords = eq.GetMesh().GetCoordinate();
  ////.GetData().AsArrayHandle<ArrayHandle<Vec3f>>();

  //vtkm::cont::Invoker{}(
  //  FiniteVolumeWorklet{ current_time }, eq.GetMesh().GetCellSet(), coords, cell_links, sln, rhs);

  //for (size_t i = 0; i < rhs.GetNumberOfValues(); i++)
  //{
  //  std::cout << "i: " << rhs.ReadPortal().Get(i) << std::endl;
  //}

  //ArrayHandle<Real> volume;
  //vtkm::cont::Invoker{}(VolumeWorklet{}, eq.GetMesh().GetCellSet(), coords, cell_links, volume);

  //std::cout << volume.ReadPortal().Get(80) << std::endl;

  //console::Error("here");
}


}
